#include "base.h"

int main(int argc,char **argv)
{




//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          //////////////////

          char dir[1024];
	  char dir_minus[1024];
	  char dir_plus[1024];


          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
	    if(strncmp(argv[i],"-dir_plus",9)==0)  sscanf(argv[i],"-dir_plus=%s",dir_plus);
            if(strncmp(argv[i],"-dir_minus",10)==0)  sscanf(argv[i],"-dir_minus=%s",dir_minus);
          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);

  // IO

   char file[1024];
   char dir_qav[1024];
   sprintf(dir_qav,"%s/QAV/",dir);

    //////////////////////////////////
    // Read Zq 
     /////////////////////////////

    double ***Zq=malloc(sizeof(double**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) Zq[imom]=malloc(sizeof(double*)*(3));
    for (int imom=0; imom<Nmom_eq; imom++) for(int ipm=0; ipm<3; ipm++) Zq[imom][ipm]=malloc(sizeof(double)*(Nclust+1));

    O55 ***DpPC=malloc(sizeof(O55**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) DpPC[imom]=malloc(sizeof(O55*)*(3));
    for (int imom=0; imom<Nmom_eq; imom++) for(int ipm=0; ipm<3; ipm++) DpPC[imom][ipm]=malloc(sizeof(O55)*(Nclust+1));

    double *aux_clust=malloc(sizeof(double)*(Nclust+1));

    for (int ipm=0; ipm<2; ipm++)
    {
   
    if (ipm==0) sprintf(dir_qav,"%s/QAV/",dir_plus);
    if (ipm==1) sprintf(dir_qav,"%s/QAV/",dir_minus);
 
    sprintf(file,"Zq");
    for (int imom=0; imom<Nmom_eq; imom++) {
	    read_input_onemomenta(dir_qav,file,aux_clust,imom);
	    for (int iclust=0; iclust<Nclust+1; iclust++) Zq[imom][ipm][iclust]=aux_clust[iclust];
     }

    for (int imom=0; imom<Nmom_eq; imom++) {
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
          {
		sprintf(file,"DpPC%d%d",iop1+1,iop2+1);
		read_input_onemomenta(dir_qav,file,aux_clust,imom);
		for (int iclust=0; iclust<Nclust+1; iclust++) DpPC[imom][ipm][iclust][iop1][iop2]=aux_clust[iclust];
	  }
    }    
    }


   //////////////////////
   //PROMEDIUM
   ///////////////////
   for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) 
   {
	Zq[imom][2][iclust]=(Zq[imom][0][iclust]+Zq[imom][1][iclust])/2.;
	for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) DpPC[imom][2][iclust][iop1][iop2]=(DpPC[imom][0][iclust][iop1][iop2]+DpPC[imom][1][iclust][iop1][iop2])/2.;
   }
     ///////////////////////
     /// WRITING OUTPUT
     ///////////////////////

      //printf(" WRITING OUTPUT IN %s\n",dir);
     sprintf(dir_qav,"%s/QAV/",dir);
     double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));

     for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq[imom][2][iclust];
     sprintf(file,"Zq");
     write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);

     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=DpPC[imom][2][iclust][iop1][iop2];
         sprintf(file,"DpPC%d%d",iop1+1,iop2+1);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
    }


    free(Zq); free(DpPC);
    return 0;

}

